## Analysis of the experimental data for Levendusky & Stecula (Cambridge Element)
## Original: Fall 2019
## This version: Winter 2021  

## This file replicates Tables 3 - 7 in the paper, plus Tables A1 - A8 in the appendix 
## Tables 1 and 2 come from the NES data, see levendusky_stecula_nes_replication.r 
## Table A9 in the appendix comes from the Lucid data, see levendusky_stecula_article_only_replication.r 

## Required R libraries 
library(tidyverse) 
library(janitor)
library(stringr)
library(lubridate)
library(lmtest)
library(plm)
library(car) 
library(stargazer)
library(modelsummary)
library(estimatr)

## Set working directory & load in the data 
# setwd("SET YOUR WORKING DIRECTORY HERE")
group <- read_csv(file="anonymized_data.csv") %>% 
  clean_names(.) 

###################################
## Data Preparation, Pre-Models  ##
################################### 

## Treatment assignment 
## 0 is control 
## 1 is polarization 
## 2 is moderation 
group <- group %>% 
  mutate(treatment = ifelse(article == "Beach",0,
                            ifelse(article == "Polarization",1,
                                   ifelse(article == "Consensus" | article == "Concensus",2,NA)))) 
## dummies for treatment 
group <- group %>% 
  mutate(control = ifelse(treatment == 0,1,0),
         polarization = ifelse(treatment == 1,1,0),
         consensus = ifelse(treatment == 2,1,0))

## For standard errors: cluster by group 
group <- group %>% 
  mutate(cluster_group = paste(session_number,group_number,sep="-"))

## Group size 
table(group$cluster_group[group$treatment==2]) ## cross-party groups 
table(group$cluster_group[group$treatment !=2])

## variable to give you the number of people in each group 
group <- group %>% 
  group_by(cluster_group) %>% 
  mutate(group_size = n())

## dummy for party ID 
group <- group %>% 
  mutate(dem = ifelse(party_id < 4,1,0))

## sample demographics 
table(group$educ) 
table(group$income) 
table(group$race)
table(group$age)
table(group$pol_interest)

## dummies for race 
group <- group %>% 
  mutate(white = ifelse(race == 1,1,0),
         black = ifelse(race == 2,1,0),
         asian = ifelse(race == 3,1,0),
         hispanic = ifelse(race == 4,1,0),
         native = ifelse(race == 5,1,0),
         other = ifelse(race == 6,1,0))

###########################################################################
## Table 3: Does Heterogeneous Discussion Reduce Affective Polarization? ##
###########################################################################

## All variables are scaled on the [0,1] interval for ease of comparison 
## As well as making a scale 

## effect on feeling thermometer ratings?  
group <- group %>% 
  mutate(other_ft = ifelse(party_id < 4,rep_ft,
                           ifelse(party_id >4,dem_ft, NA)),
         same_ft = ifelse(party_id < 4, dem_ft,
                          ifelse(party_id > 4,rep_ft, NA)),
         ft_diff = same_ft - other_ft, 
         other_cand = ifelse(party_id < 4, djtft,
                             ifelse(party_id > 4, hrcft, NA)),
         same_cand = ifelse(party_id < 4, hrcft,
                             ifelse(party_id > 4, djtft, NA)), 
         other_ft = other_ft/100,
         same_ft = same_ft/100,
         ft_diff = ft_diff/100,
         other_cand = other_cand/100,
         same_cand = same_cand/100, 
         cand_diff = same_cand - other_cand) 

m.oft <- lm_robust(other_ft ~ as.factor(treatment) + as.factor(session_number),
                data = group,
                cluster = cluster_group)
m.ocft <- lm_robust(other_cand ~ as.factor(treatment) + as.factor(session_number),
                    data = group,
                    cluster = cluster_group)

## ft modifications 
## 10 degrees is the 50th percentile in the control condition 
group <- group %>% 
  mutate(less10 = ifelse(other_ft  < .11,1,0),
         gtfifty = ifelse(other_ft > .49,1,0))

m.oftlt <- lm_robust(less10 ~ as.factor(treatment) + as.factor(session_number),
                data = group,
                cluster = cluster_group)
m.oftgtf <- lm_robust(gtfifty ~ as.factor(treatment) + as.factor(session_number),
                data = group,
                cluster = cluster_group)
## one striking difference: heterog. vs. control gives a +10 degree FT diff for the other party (20 degrees vs. 30 degrees, 50\% increase)

## Trait Ratings 
## combined trait ratings 
rev_hyp <- -1*group$other_hypocritical + 6 
rev_self <- -1*group$other_selfish + 6 
rev_mean <- -1*group$other_mean + 6 
combo_traits <- cbind(group$other_american,
                      group$other_intelligent,
                      group$other_honest,
                      group$other_generous,
                      group$other_open,
                      rev_hyp,
                      rev_self,
                      rev_mean) 
group$avg_trait <- apply(combo_traits,1,function(x){mean(x,na.rm=T)})
group$avg_trait <- group$avg_trait/5 
m.tr <- lm_robust(avg_trait ~ as.factor(treatment) + as.factor(session_number),
                data = group, 
                subset = party_id != 4,
                cluster = cluster_group)

## alpha of combined trait rating scale? 
foralpha <- cbind(group$other_american,
                  group$other_intelligent,
                  group$other_honest,
                  group$other_generous,
                  group$other_open,
                  group$other_selfish,
                  group$other_hypocritical,
                  group$other_mean)
psych::alpha(foralpha, check.keys = T) 
## alpha = 0.83 

## Trust in the other party 
group <- group %>% 
  mutate(trust_other = ifelse(party_id < 4,trust_rep,
                ifelse(party_id >4,trust_dem, NA)),
         trust_same = ifelse(party_id < 4, trust_dem,
                          ifelse(party_id > 4,trust_rep, NA)),
         trust_other = trust_other/5,
         trust_same = trust_same/5, 
         trust_diff = trust_same - trust_other) 

m.trust <- lm_robust(trust_other ~ as.factor(treatment) + as.factor(session_number),
                 data = group,
                 cluster = cluster_group)

## social distance measures 
group <- group %>% 
  mutate(rev_friend = (-1*comfort_friend) + 5,
         rev_neighbor = (-1*comfort_neighbor)+5,
         rev_discuss = (-1*comfort_discuss) + 5) %>% 
  rowwise(.) %>% 
  mutate(soc_dist = mean(c(rev_friend,
                           rev_neighbor,
                           rev_discuss,
                           comfort_marriage), na.rm=T)) %>% 
  ungroup(.) 
group$soc_dist <- group$soc_dist/4 

## do these scale? 
psych::alpha(cbind(group$comfort_marriage, 
                   group$rev_discuss,
                   group$rev_friend,
                   group$rev_neighbor))
## alpha = 0.8, one scale 

m.socdist <- lm_robust(soc_dist ~ as.factor(treatment) + as.factor(session_number),
                data = group, 
                subset = (party_id != 4),
                cluster = cluster_group)

## make into an overall index 
## other party FT, other-party trust, other_party traits, social distance 
psych::alpha(cbind(group$other_ft,
                   group$avg_trait,
                   group$trust_other,
                   group$soc_dist)) 
## alpha = 0.75 

group <- group %>% 
  rowwise(.) %>% 
  mutate(out_party_affect = mean(c(other_ft,avg_trait,trust_other,soc_dist),na.rm=T)) %>% 
  ungroup(.) 

m.index <- lm_robust(out_party_affect ~ as.factor(treatment) + as.factor(session_number),
                     data = group, 
                     subset = (party_id != 4),
                     cluster = cluster_group) 


## output as a table 
## need to do the goofy fix noted here: https://declaredesign.org/r/estimatr/articles/regression-tables.html
## make non-robust regs, then pass to starprep function 
mn.index <-lm(out_party_affect ~ as.factor(treatment) + as.factor(session_number),
              data = group,
              subset = party_id != 4)
mn.oft <-lm(other_ft ~ as.factor(treatment) + as.factor(session_number),
            data = group)
mn.ocft <- lm(other_cand ~ as.factor(treatment) + as.factor(session_number),
              data = group)
mn.oftlt <- lm(less10 ~ as.factor(treatment) + as.factor(session_number),
              data = group)
mn.oftgtf <- lm(gtfifty ~ as.factor(treatment) + as.factor(session_number),
               data = group)
mn.tr <-lm(avg_trait ~ as.factor(treatment) + as.factor(session_number),
            data = group,
            subset = party_id != 4)
mn.trust <- lm(trust_other ~ as.factor(treatment) + as.factor(session_number),
               data = group) 
mn.socdist <- lm(soc_dist ~ as.factor(treatment) + as.factor(session_number),
                 subset = (party_id != 4),
                 data = group)

stargazer(mn.index,mn.oft,mn.oftlt,mn.oftgtf,mn.ocft,mn.tr,mn.trust,mn.socdist,
          se = starprep(m.index,m.oft,m.oftlt,m.oftgtf,m.ocft,m.tr,m.trust,m.socdist),
          column.labels = c("Index","Out-Party FT","Out-Party FT < 10",
                            "Out-Party FT >50","Out-Party Cand. FT",
                            "Out-Party Traits","Out-Party Trust","Social Distance"),
          covariate.labels = c("Homogeneous Discussion","Heterogeneous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/table3_main_effects.html")

##########################################################
## Table 4: Heterogeneous Effects by Partisan Strength? ##
##########################################################

## Effects separated by partisan strength? 
group <- group %>% 
  mutate(sp = ifelse(party_id == 1 | party_id == 7,1,0), 
         pid_strength = ifelse(party_id == 4,0,
                               ifelse(party_id == 3 | party_id == 5,1,
                                      ifelse(party_id == 2 | party_id == 6,2, 
                                             ifelse(party_id == 1 | party_id == 7,3,NA)))),
         sorted = ifelse((party_id < 4 & pol_views < 4) | (party_id > 4 & pol_views > 4),1,0))

m16 <- lm_robust(out_party_affect ~ as.factor(treatment)*sp + as.factor(session_number),
                 data = group,
                 subset = party_id != 4,
                 cluster = cluster_group) 
lht(m16, c("as.factor(treatment)2 + as.factor(treatment)2:sp = 0"))
m17 <- lm_robust(other_ft ~ as.factor(treatment)*sp + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
lht(m17, c("as.factor(treatment)2 + as.factor(treatment)2:sp = 0"))
m19 <- lm_robust(other_cand ~ as.factor(treatment)*sp + as.factor(session_number),
                 data = group, 
                 cluster = cluster_group) 
lht(m19, c("as.factor(treatment)2 + as.factor(treatment)2:sp = 0"))
m20 <- lm_robust(trust_other ~ as.factor(treatment)*sp + as.factor(session_number),
                 data = group,
                 cluster = cluster_group)
lht(m20, c("as.factor(treatment)2 + as.factor(treatment)2:sp = 0"))
m21 <- lm_robust(avg_trait ~ as.factor(treatment)*sp + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group)
lht(m21, c("as.factor(treatment)2 + as.factor(treatment)2:sp = 0"))
m22 <- lm_robust(soc_dist ~ as.factor(treatment)*sp + as.factor(session_number),
                 data = group, 
                 subset = party_id != 4,
                 cluster = cluster_group)
lht(m22, c("as.factor(treatment)2 + as.factor(treatment)2:sp = 0"))

## output to a table 
models <- list()
models[["AP Index"]] <- m16  
models[["Out-Party FT"]] <- m17 
models[["Out-Party Cand. FT"]] <- m19 
models[["Out-Party Traits"]] <- m21
models[["Out-Party Trust"]] <- m20 
models[["Social Distance (Reversed)"]] <- m22 

msummary(models,
         coef_map = c("(Intercept)" = "Constant",
                      "as.factor(treatment)1" = "Homogeneous Discussion",
                      "as.factor(treatment)2" = "Heterogeneous Discussion",
                      "sp" = "Strong Partisan",
                      "as.factor(treatment)1:sp" = "Strong Partisan*Homog. Discussion",
                      "as.factor(treatment)2:sp" = "Strong Partisan*Heterog. Discussion"),
         stars = T, 
         output =  "tables/table4_het_efx_pstrength.html")


###################################################################
## Table 5: Heterogeneous Effects by Partisan Media Consumption? ##
###################################################################

## For Republicans: avg # of days of Fox & right-wing websites 
## For Democrats: avg # of days of MSNBC & left-wing websites 

group <- group %>% 
  rowwise(.) %>% 
  mutate(lw_media = mean(c(msnbc,left_lean_online)),
         rw_media = mean(c(fox_news,right_lean_online))) %>% 
  ungroup(.) 

group <- group %>% 
  mutate(partisan_media = ifelse(party_id < 4,lw_media,
                                 ifelse(party_id > 4, rw_media,NA)))

m23 <- lm_robust(out_party_affect ~ as.factor(treatment)*partisan_media + as.factor(session_number),
                 data = group,
                 subset = party_id != 4,
                 cluster = cluster_group) 
lht(m23, c("as.factor(treatment)2 + as.factor(treatment)2:partisan_media = 0"))
m24 <- lm_robust(other_ft ~ as.factor(treatment)*partisan_media + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
lht(m24, c("as.factor(treatment)2 + as.factor(treatment)2:partisan_media = 0"))
m25 <- lm_robust(other_cand ~ as.factor(treatment)*partisan_media + as.factor(session_number),
                 data = group, 
                 cluster = cluster_group) 
lht(m25, c("as.factor(treatment)2 + as.factor(treatment)2:partisan_media = 0"))
m26 <- lm_robust(trust_other ~ as.factor(treatment)*partisan_media + as.factor(session_number),
                 data = group,
                 cluster = cluster_group)
lht(m26, c("as.factor(treatment)2 + as.factor(treatment)2:partisan_media = 0"))
m27 <- lm_robust(avg_trait ~ as.factor(treatment)*partisan_media + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group)
lht(m27, c("as.factor(treatment)2 + as.factor(treatment)2:partisan_media = 0"))
m28 <- lm_robust(soc_dist ~ as.factor(treatment)*partisan_media + as.factor(session_number),
                 data = group, 
                 subset = party_id != 4,
                 cluster = cluster_group)
lht(m28, c("as.factor(treatment)2 + as.factor(treatment)2:partisan_media = 0"))

## output to a table 
models <- list()
models[["AP Index"]] <- m23  
models[["Out-Party FT"]] <- m24 
models[["Out-Party Cand. FT"]] <- m25 
models[["Out-Party Traits"]] <- m27
models[["Out-Party Trust"]] <- m26 
models[["Social Distance (Reversed)"]] <- m28 

msummary(models,
         coef_map = c("(Intercept)" = "Constant",
                      "as.factor(treatment)1" = "Homogeneous Discussion",
                      "as.factor(treatment)2" = "Heterogeneous Discussion",
                      "partisan_media" = "Partisan Media Consumption",
                      "as.factor(treatment)1:partisan_media" = "Media*Homog. Discussion",
                      "as.factor(treatment)2:partisan_media" = "Media*Heterog. Discussion"),
         stars = T, 
         output =  "tables/table5_het_efx_part_media.html")
 

##########################
## Table 6: Mechanisms? ## 
##########################

psych::alpha(cbind(group$rev_discuss,
                   group$understand_other,
                   group$reasonable_position,
                   group$respect_opinion))
## these do not scale! Analyze them separately 

mech1 <- lm_robust(understand_other ~ as.factor(treatment) + as.factor(session_number),
                   subset = party_id != 4, 
                   data = group, 
                   cluster = cluster_group)
mech2 <- lm_robust(fact_evidence ~ as.factor(treatment) + as.factor(session_number),
                   subset = party_id != 4, 
                   data = group,
                   cluster = cluster_group) 
mech3 <- lm_robust(reasonable_position ~ as.factor(treatment) + as.factor(session_number), 
                   data = group,
                   subset = party_id != 4, 
                   cluster = cluster_group) 
mech4 <- lm_robust(respect_opinion ~ as.factor(treatment) + as.factor(session_number), 
                   data = group,
                   subset = party_id != 4, 
                   cluster = cluster_group)
mechn1 <- lm(understand_other ~ as.factor(treatment) + as.factor(session_number),
                   subset = party_id != 4, 
                   data = group)
mechn2 <- lm(fact_evidence ~ as.factor(treatment) + as.factor(session_number),
                   subset = party_id != 4, 
                   data = group) 
mechn3 <- lm(reasonable_position ~ as.factor(treatment) + as.factor(session_number), 
                   data = group,
                   subset = party_id != 4) 
mechn4 <- lm(respect_opinion ~ as.factor(treatment) + as.factor(session_number), 
                   data = group,
                   subset = party_id != 4)


## Perceptions of Common Ground? 
cg <- cbind(group$common_ground, 
            group$agree_more)
psych::alpha(cg)
## alpha = 0.72 
group$avg_cg <- apply(cg,1,function(x){mean(x,na.rm=T)})
group$avg_cg <- group$avg_cg/5

mech.cg <- lm_robust(avg_cg ~ as.factor(treatment) + as.factor(session_number),
                      data = group, 
                     subset = party_id != 4, 
                     cluster = cluster_group) 
mn.cg <- lm(avg_cg ~ as.factor(treatment) + as.factor(session_number),
           data = group, 
           subset = party_id != 4) 

stargazer(mn.cg,mechn1,mechn3,mechn4,
          se = starprep(mech.cg,mech1,mech3,mech4),
          column.labels = c("Common Ground","Understand Other","Reasonable Beliefs","Other Respects You"),
          covariate.labels = c("Homogeneous Discussion","Heterogeneous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/table6_mechanisms.html")

################################################################
## Table 7: Do the Effects Endure? (Follow-Up Study Analysis) ##
################################################################

## Who takes the follow-up study? 
group <- group %>% 
  mutate(take_followup = ifelse(!is.na(q3),1,0),
         white = ifelse(race == 1,1,0), 
         dem = ifelse(party_id < 4,1,0))

mean(group$take_followup,na.rm=T)


## Effects on FT Ratings, Trust, PID Importance, Common Ground   
group <- group %>% 
  mutate(follow_other_ft = ifelse(party_id < 4, q7_2,
                                  ifelse(party_id > 4, q7_1,NA)),
         follow_same_ft = ifelse(party_id > 4, q7_2,
                                 ifelse(party_id < 4, q7_1,NA)),
         follow_ft_diff = follow_same_ft - follow_other_ft,
         follow_other_cand = ifelse(party_id < 4, q7_4,
                                  ifelse(party_id > 4, q7_3,NA)),
         follow_same_cand = ifelse(party_id > 4, q7_4,
                                 ifelse(party_id < 4, q7_3,NA)),
         follow_other_trust = ifelse(party_id > 4, q26,
                                     ifelse(party_id < 4, q27,NA)),
         follow_pid_import = ifelse(party_id > 4, q5,
                                    ifelse(party_id < 4, q4,NA)))
## rescale to [0,1] 
group$follow_other_ft <- group$follow_other_ft/100 
group$follow_other_cand <- group$follow_other_cand/100 
group$follow_other_trust <- group$follow_other_trust/5 

## make the follow-up common ground measure 
group <- group %>% 
  rowwise(.) %>% 
  mutate(follow_cg = mean(c(q16,q17),na.rm=T)) 

## combined trait ratings (follow-up)
group <- group %>% 
  mutate(follow_american = ifelse(!is.na(q28_1),(-1*q28_1)+6,
                                  ifelse(!is.na(q29_1),(-1*q29_1)+6,NA)),
         follow_intelligent = ifelse(!is.na(q28_2),(-1*q28_2)+6,
                                     ifelse(!is.na(q29_2),(-1*q29_2)+6,NA)),
         follow_honest = ifelse(!is.na(q28_3),(-1*q28_3)+6,
                                     ifelse(!is.na(q29_3),(-1*q29_3)+6,NA)),
         follow_open = ifelse(!is.na(q28_4),(-1*q28_4)+6,
                                     ifelse(!is.na(q29_4),(-1*q29_4)+6,NA)),
         follow_generous = ifelse(!is.na(q28_5),(-1*q28_5)+6,
                                     ifelse(!is.na(q29_5),(-1*q29_5)+6,NA)),
         follow_hypocritical = ifelse(!is.na(q28_6),q28_6,
                                       ifelse(!is.na(q29_6),q29_6,NA)),
         follow_selfish = ifelse(!is.na(q28_7),q28_7,
                                      ifelse(!is.na(q29_7),q29_6,NA)),
         follow_mean = ifelse(!is.na(q28_8),q28_8,
                                      ifelse(!is.na(q29_8),q29_6,NA))) 

follow_combo_traits <- cbind(group$follow_american,
                             group$follow_intelligent,
                             group$follow_honest,
                             group$follow_generous,
                             group$follow_open,
                             group$follow_hypocritical,
                             group$follow_selfish,
                             group$follow_mean)
group$follow_avg_trait <- apply(follow_combo_traits,1,function(x){mean(x,na.rm=T)})
group$follow_avg_trait <- group$follow_avg_trait/5 

## Social Distance Measures 
group <- group %>% 
  mutate(follow_comfort_discuss = ifelse(!is.na(q18),(-1*q18)+5,
                                           ifelse(!is.na(q19),(-1*q19)+5,NA)),
         follow_comfort_neighbor = ifelse(!is.na(q20),(-1*q20)+5,
                                           ifelse(!is.na(q21),(-1*q21)+5,NA)),
         follow_comfort_friends = ifelse(!is.na(q22),(-1*q22)+5,
                                            ifelse(!is.na(q23),(-1*q23)+5,NA)),
         follow_comfort_marriage = ifelse(!is.na(q24),(-1*q24)+5,
                                           ifelse(!is.na(q25),(-1*q25)+5,NA))) %>% 
  rowwise(.) %>% 
  mutate(follow_social_dist = mean(c(follow_comfort_discuss,
                                     follow_comfort_neighbor,
                                     follow_comfort_friends,
                                     follow_comfort_marriage), na.rm = T)) %>% 
  ungroup(.) 
group$follow_social_dist <- group$follow_social_dist/4 

## Index of all measures 
group <- group %>% 
  rowwise(.) %>% 
  mutate(follow_out_party_affect = mean(c(follow_other_ft,
                                          follow_avg_trait,
                                          follow_other_trust,
                                          follow_social_dist), na.rm=T)) %>%
  ungroup(.) 
           

## Regression Models 
mf.all <- lm_robust(follow_out_party_affect ~ as.factor(treatment) + as.factor(session_number),
                    data = group, 
                    subset = party_id != 4, 
                    cluster = cluster_group)  
mf.ft <- lm_robust(follow_other_ft ~ as.factor(treatment) + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
mf.ftcand <- lm_robust(follow_other_cand ~ as.factor(treatment) + as.factor(session_number),
                 data = group, 
                 cluster = cluster_group) 
mf.trust <- lm_robust(follow_other_trust ~ as.factor(treatment) + as.factor(session_number),
                 data = group, 
                 cluster = cluster_group) 
mf.trait <- lm_robust(follow_avg_trait ~ as.factor(treatment) + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group) 
mf.socdist <- lm_robust(follow_social_dist ~ as.factor(treatment) + as.factor(session_number),
                        subset = party_id != 4, 
                         data = group, 
                         cluster = cluster_group) 

## non-robust version for stargazer 
mnf.all <- lm(follow_out_party_affect ~ as.factor(treatment) + as.factor(session_number),
              subset = party_id != 4, 
              data = group)  
mnf.ft <- lm(follow_other_ft ~ as.factor(treatment) + as.factor(session_number),
                 data = group) 
mnf.ftcand <- lm(follow_other_cand ~ as.factor(treatment) + as.factor(session_number),
                 data = group) 
mnf.trust <- lm(follow_other_trust ~ as.factor(treatment) + as.factor(session_number),
                 data = group) 
mnf.trait <- lm(follow_avg_trait ~ as.factor(treatment) + as.factor(session_number),
                 data = group,
                 subset = party_id != 4) 
mnf.socdist <- lm(follow_social_dist ~ as.factor(treatment) + as.factor(session_number),
                  subset = party_id != 4, 
                 data = group) 


## output to a table 
stargazer(mnf.all,mnf.ft,mnf.ftcand,mnf.trait,mnf.trust,mnf.socdist,
          se = starprep(mf.all,mf.ft,mf.ftcand,mf.trait,mf.trust,mf.socdist),
          column.labels = c("AP Index","Out-Party FT","Out-Party Cand. FT",
                            "Out-Party Traits","Out-Party Trust","Social Distance (Reversed)"),
          covariate.labels = c("Homogeneous Discussion","Heterogenenous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/table7_duration_effects.html")

##################################
## Table 8: Downstream Effects? ##
##################################

## Future Group 
group <- group %>% 
  mutate(homog_group = ifelse(party_id < 4, group_type,
                              ifelse(party_id > 4, (-1*group_type)+8,NA)),
         dem = ifelse(party_id < 4,1,0))

m14 <- lm_robust(homog_group ~ as.factor(treatment) + as.factor(session_number),
                 data = group, 
                 subset = party_id != 4, 
                 cluster = cluster_group) 
## no effect overall 
m15 <- lm_robust(group_type ~ as.factor(treatment)*sp + as.factor(session_number),
          data = group,
          subset = party_id != 4, 
          cluster = cluster_group)  
lht(m15, c("as.factor(treatment)2 + as.factor(treatment)2:sp = 0"))
## effect for weak/leaning partisans (p < .1) --> they prefer a more heterogeneous group 

## Partisan Identity importance 
m16 <- lm_robust(pid_import ~ as.factor(treatment) + as.factor(session_number),
                  data = group,
                  subset = party_id != 4, 
                  cluster = cluster_group)  
m17 <- lm_robust(pid_import ~ as.factor(treatment)*sp + as.factor(session_number),
                  data = group,
                  subset = party_id != 4, 
                  cluster = cluster_group)

## Output to a Table 
models <- list()
models[["Future Discussion (1)"]] <- m14 
models[["PID Importance (1)"]] <- m16  
models[["Future Discussion (2)"]] <- m15
models[["PID Importance (2)"]] <- m17 


msummary(models,
         coef_map = c("(Intercept)" = "Constant",
                      "as.factor(treatment)1" = "Homogeneous Discussion",
                      "as.factor(treatment)2" = "Heterogeneous Discussion",
                      "sp" = "Strong Partisan",
                      "as.factor(treatment)1:sp" = "Strong Partisan*Homog. Discussion",
                      "as.factor(treatment)2:sp" = "Strong Partisan*Heterog. Discussion"),
         stars = T, 
         output =  "tables/table8_downstream_efx.html")


#############################
## Tables for the Appendix ## 
#############################

## Table A1: heterogeneity based on group size? 
ms0 <- lm_robust(out_party_affect ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group) 
ms1 <- lm_robust(other_ft ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
ms2 <- lm_robust(same_ft ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
ms3 <- lm_robust(ft_diff ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
ms4 <- lm_robust(less10 ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
ms5 <- lm_robust(gtfifty ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
ms6 <- lm_robust(other_cand ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
ms7 <- lm_robust(trust_other ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
ms8 <- lm_robust(avg_trait ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group) 
ms9 <- lm_robust(soc_dist ~ as.factor(treatment)*group_size + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group) 

## non-robust versions for stargazer 
mns0 <- lm(out_party_affect ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group, 
           subset = party_id != 4) 
mns1 <- lm(other_ft ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group) 
mns2 <- lm(same_ft ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group) 
mns3 <- lm(ft_diff ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group) 
mns4 <- lm(less10 ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group) 
mns5 <- lm(gtfifty ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group) 
mns6 <- lm(other_cand ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group) 
mns7 <- lm(trust_other ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group) 
mns8 <- lm(avg_trait ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group,
           subset = party_id != 4)
mns9 <- lm(soc_dist ~ as.factor(treatment)*group_size + as.factor(session_number),
           data = group,
           subset = party_id != 4)

## output as an appendix table 
stargazer(mns0,mns1,mns4,mns5,mns6,mns8,mns7,mns9,
          se = starprep(ms0,ms1,ms4,ms5,ms6,ms7,ms8,ms9),
          column.labels = c("AP Index","Out-Party FT","Out-Party FT < 10",
                            "Out-Party FT >50","Out-Party Cand. FT",
                            "Out-Party Traits","Out-Party Trust",
                            "Social Distance (Reversed)"),
          covariate.labels = c("Homogeneous Discussion","Heterogeneous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/appendix_het_efx_size.html")

## Table A2: with the same-party/differenced feeling thermometer ratings 
m.sft <- lm_robust(same_ft ~ as.factor(treatment) + as.factor(session_number),
                   data = group,
                   cluster = cluster_group)
m.ftd <- lm_robust(ft_diff ~ as.factor(treatment) + as.factor(session_number),
                   data = group,
                   cluster = cluster_group)
m.scft <- lm_robust(same_cand ~ as.factor(treatment) + as.factor(session_number),
                    data = group,
                    cluster = cluster_group) 
m.cftd <- lm_robust(cand_diff ~ as.factor(treatment) + as.factor(session_number),
                    data = group,
                    cluster = cluster_group) 
m.st <- lm_robust(trust_same ~ as.factor(treatment) + as.factor(session_number),
                  data = group,
                  cluster = cluster_group)
m.td <- lm_robust(trust_diff ~ as.factor(treatment) + as.factor(session_number),
                  data = group,
                  cluster = cluster_group) 

## non-robust version for stargazer 
mn.sft <-lm(same_ft ~ as.factor(treatment) + as.factor(session_number),
            data = group)
mn.ftd <-lm(ft_diff ~ as.factor(treatment) + as.factor(session_number),
            data = group)
mn.scft <- lm(same_cand ~ as.factor(treatment) + as.factor(session_number),
              data = group)
mn.cftd <- lm(cand_diff ~ as.factor(treatment) + as.factor(session_number),
              data = group) 
mn.st <- lm(trust_same ~ as.factor(treatment) + as.factor(session_number),
            data = group)
mn.td <- lm(trust_diff ~ as.factor(treatment) + as.factor(session_number),
            data = group)

stargazer(mn.ftd,mn.cftd,mn.td,
          se = starprep(m.sft,m.ftd,m.scft,m.cftd,m.st,m.td),
          column.labels = c("Party FT Difference",
                            "Cand. FT Difference",
                            "Difference in Trust"),
          covariate.labels = c("Homogeneous Discussion","Heterogeneous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/appendix_ft_diff_measures.html") 

## Table A3: positive/negative trait ratings 
positive_traits <- cbind(group$other_american,
                         group$other_intelligent,
                         group$other_honest,
                         group$other_open,
                         group$other_generous)
group$avg_positive <- apply(positive_traits,1,function(x){mean(x,na.rm=T)})
library(psych) 
psych::alpha(positive_traits)
m.tp <- lm_robust(avg_positive ~ as.factor(group$treatment) + as.factor(session_number),
                  data = group,
                  subset = party_id != 4, 
                  cluster = cluster_group) 
mn.tp <- lm(avg_positive ~ as.factor(group$treatment) + as.factor(session_number),
            data = group,
            subset = party_id != 4) 
## find the expected effects 

## negative traits 
negative_traits <- cbind(group$other_hypocritical,
                         group$other_selfish,
                         group$other_mean) 
group$avg_neg <- apply(negative_traits,1,function(x){mean(x,na.rm=T)})
psych::alpha(negative_traits) 
m.tn <- lm_robust(avg_neg ~ as.factor(treatment) + as.factor(session_number),
                  data = group,
                  subset = party_id != 4, 
                  cluster = cluster_group)
mn.tn <- lm(avg_neg ~ as.factor(treatment) + as.factor(session_number),
            data = group,
            subset = party_id != 4) 
## effect is not sig, consistent w/ Mullinix arg on negative >> positive 

stargazer(mn.tp,mn.tn,
          se = starprep(m.tp,m.tn),
          column.labels = c("Positive Traits","Negative Traits"),
          covariate.labels = c("Homogeneous Discussion","Heterogeneous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/appendix_pos_neg_traits.html") 

## Table A4: individual social distance items 
m.cm <- lm_robust(comfort_marriage ~ as.factor(treatment) + as.factor(session_number),
                  data = group, 
                  subset = (  party_id != 4),
                  cluster = cluster_group)
m.cd <- lm_robust(rev_discuss ~ as.factor(treatment) + as.factor(session_number),
                  data = group, 
                  subset = (  party_id != 4),
                  cluster = cluster_group)
m.cf <- lm_robust(rev_friend ~ as.factor(treatment) + as.factor(session_number),
                  data = group, 
                  subset = (  party_id != 4),
                  cluster = cluster_group)
m.cn <- lm_robust(rev_neighbor ~ as.factor(treatment) + as.factor(session_number),
                  data = group, 
                  subset = (  party_id != 4),
                  cluster = cluster_group)
mn.cm <- lm(comfort_marriage ~ as.factor(treatment) + as.factor(session_number),
            data = group, 
            subset = (  party_id != 4))
mn.cd <- lm(rev_discuss ~ as.factor(treatment) + as.factor(session_number),
            data = group, 
            subset = (  party_id != 4))
mn.cf <- lm(rev_friend ~ as.factor(treatment) + as.factor(session_number),
            data = group, 
            subset = (  party_id != 4))
mn.cn <- lm(rev_neighbor ~ as.factor(treatment) + as.factor(session_number),
            data = group, 
            subset = (  party_id != 4))

stargazer(mn.cm,mn.cd,mn.cf,mn.cn,
          se = starprep(m.cm,m.cd,m.cf,m.cn),
          column.labels = c("Marriage","Discussion","Friendship","Neighbors"),
          covariate.labels = c("Homogeneous Discussion","Heterogeneous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/appendix_soc_dist.html") 

## Table A5: who takes the follow-up study? 
m.pred0 <- glm(take_followup ~ polarization + consensus,
               family=binomial(link=probit),
               data=group)
m.pred1 <- glm(take_followup ~ polarization + consensus + dem + sp + pol_views + pol_interest + pol_activity + days_talk,
               family = binomial(link=probit),
               na.action = "na.exclude",
               data = group) 
m.pred2 <- glm(take_followup ~ polarization + consensus + dem + sp + pol_views + pol_interest + pol_activity + days_talk 
               + educ + income + black + asian + hispanic + native + other + age + female,  
               family = binomial(link=probit),
               na.action = "na.exclude",
               data = group) 
summary(m.pred0)
summary(m.pred1)
summary(m.pred2)
## No treatment/party/strong partisan/interest effects 
## Effects for better educated, younger, and women (those folks are more likley to take our study) 

## table of participating in our follow-up study 
stargazer(m.pred0, m.pred1, m.pred2,
          column.labels = c("Model 1", "Model 2", "Model 3"), 
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/appendix_table_follow_participation.html")

## Table A6: heterogeneous effects based on days since the original study? 
group <- group %>% 
  mutate(inperson_date = dmy(session_date), 
         tmp = mdy_hm(end_date), 
         follow_up_date = as.Date(tmp),
         days_since_study = as.numeric(follow_up_date - inperson_date),
         ## one weird recording errror, fix 
         days_since_study = ifelse(days_since_study <0,NA,days_since_study))

table(group$days_since_study)

md0 <- lm_robust(follow_out_party_affect ~ as.factor(treatment)*days_since_study + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group) 
md1 <- lm_robust(follow_other_ft ~ as.factor(treatment)*days_since_study + as.factor(session_number),
                 data = group,
                 cluster = cluster_group) 
md3 <- lm_robust(follow_other_cand ~ as.factor(treatment)*days_since_study + as.factor(session_number),
                 data = group, 
                 cluster = cluster_group) 
md4 <- lm_robust(follow_other_trust ~ as.factor(treatment)*days_since_study + as.factor(session_number),
                 data = group, 
                 cluster = cluster_group)
md5 <- lm_robust(follow_avg_trait ~ as.factor(treatment)*days_since_study + as.factor(session_number),
                 data = group, 
                 subset = party_id != 4, 
                 cluster = cluster_group)
md6 <- lm_robust(follow_social_dist ~ as.factor(treatment)*days_since_study + as.factor(session_number),
                 data = group, 
                 subset = party_id != 4, 
                 cluster = cluster_group)

## non-robst models for stargazer 
mnd0 <- lm(follow_out_party_affect ~ as.factor(treatment)*days_since_study + as.factor(session_number),
           data = group,
           subset = party_id != 4) 
mnd1 <- lm(follow_other_ft ~ as.factor(treatment)*days_since_study + as.factor(session_number),
           data = group) 
mnd3 <- lm(follow_other_cand ~ as.factor(treatment)*days_since_study + as.factor(session_number),
           data = group) 
mnd4 <- lm(follow_other_trust ~ as.factor(treatment)*days_since_study + as.factor(session_number),
           data = group)
mnd5 <- lm(follow_avg_trait ~ as.factor(treatment)*days_since_study + as.factor(session_number),
           data = group, 
           subset = party_id != 4)
mnd6 <- lm(follow_social_dist ~ as.factor(treatment)*days_since_study + as.factor(session_number),
           data = group, 
           subset = party_id != 4)

## table w/ delay for appendix 
stargazer(mnd0,mnd1,mnd3,mnd5,mnd4,mnd6,
          se = starprep(md0,md1,md3,md5,md4,md6),
          column.labels = c("AP Index","Out-Party FT","Out-Party Cand. FT",
                            "Out-Party Traits","Out-Party Trust","Social Distance (Reversed)"),
          covariate.labels = c("Homogeneous Discussion","Heterogenenous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/appendix_table_delay_appendix.html")

## Table A7: Control for age, education, and gender (sig. predictors of taking the study)
mc0 <- lm_robust(follow_out_party_affect ~ as.factor(treatment) + age + educ + female + as.factor(session_number), 
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group) 
mc1 <- lm_robust(follow_other_ft ~ as.factor(treatment) + age + educ + female + as.factor(session_number), 
                 data = group,
                 cluster = cluster_group) 
mc3 <- lm_robust(follow_other_cand ~ as.factor(treatment) + age + educ + female + as.factor(session_number),
                 data = group, 
                 cluster = cluster_group) 
mc4 <- lm_robust(follow_other_trust ~ as.factor(treatment) + age + educ + female + as.factor(session_number),
                 data = group, 
                 cluster = cluster_group) 
mc5 <- lm_robust(follow_avg_trait ~ as.factor(treatment) + age + educ + female + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group) 
mc6 <- lm_robust(follow_social_dist ~ as.factor(treatment) + age + educ + female + as.factor(session_number),
                 data = group,
                 subset = party_id != 4, 
                 cluster = cluster_group) 

## non-robust for stargazer 
mnc0 <- lm(follow_out_party_affect ~ as.factor(treatment) + age + educ + female + as.factor(session_number), 
           data = group,
           subset = party_id != 4) 
mnc1 <- lm(follow_other_ft ~ as.factor(treatment) + age + educ + female + as.factor(session_number), 
           data = group) 
mnc3 <- lm(follow_other_cand ~ as.factor(treatment) + age + educ + female + as.factor(session_number),
           data = group) 
mnc4 <- lm(follow_other_trust ~ as.factor(treatment) + age + educ + female + as.factor(session_number),
           data = group) 
mnc5 <- lm(follow_avg_trait ~ as.factor(treatment) + age + educ + female + as.factor(session_number),
           data = group,
           subset = party_id != 4) 
mnc6 <- lm(follow_social_dist ~ as.factor(treatment) + age + educ + female + as.factor(session_number),
           data = group,
           subset = party_id != 4)

## table w/ controls for appendix 
stargazer(mnc0, mnc1, mnc3, mnc5, mnc4, mnc6,
          se = starprep(mc0,mc1,mc3,mc5,mc4,mc6), 
          column.labels = c("AP Index","Out-Party FT","Out-Party Cand. FT",
                            "Out-Party Traits","Out-Party Trust","Social Distance (Reversed)"),
          covariate.labels = c("Homogeneous Discussion","Heterogenenous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/table_control_appendix.html")

## Table A8: Estimating Table 3 on just those who take the follow-up study 
mf.index <- lm_robust(out_party_affect ~ as.factor(treatment) + as.factor(session_number),
                      data = group, 
                      subset = (party_id != 4) & take_followup == 1, 
                      cluster = cluster_group) 
mf.oft <- lm_robust(other_ft ~ as.factor(treatment) + as.factor(session_number),
                    data = group, 
                    subset = (party_id != 4) & take_followup == 1, 
                    cluster = cluster_group) 
mf.lt <- lm_robust(less10 ~ as.factor(treatment) + as.factor(session_number),
                   data = group, 
                   subset = (party_id != 4) & take_followup == 1, 
                   cluster = cluster_group) 
mf.gtf <- lm_robust(gtfifty ~ as.factor(treatment) + as.factor(session_number),
                    data = group, 
                    subset = (party_id != 4) & take_followup == 1, 
                    cluster = cluster_group) 
mf.ocft <- lm_robust(other_cand ~ as.factor(treatment) + as.factor(session_number),
                     data = group, 
                     subset = (party_id != 4) & take_followup == 1, 
                     cluster = cluster_group) 
mf.trait <- lm_robust(avg_trait ~ as.factor(treatment) + as.factor(session_number),
                      data = group, 
                      subset = (party_id != 4) & take_followup == 1, 
                      cluster = cluster_group) 
mf.trust <- lm_robust(trust_other ~ as.factor(treatment) + as.factor(session_number),
                      data = group, 
                      subset = (party_id != 4) & take_followup == 1, 
                      cluster = cluster_group) 
mf.sd <- lm_robust(soc_dist ~ as.factor(treatment) + as.factor(session_number),
                   data = group, 
                   subset = (party_id != 4) & take_followup == 1, 
                   cluster = cluster_group) 
## results replicate, indeed get slightly larger 

## non-robust versions for stargazer 
mfn.index <- lm(out_party_affect ~ as.factor(treatment) + as.factor(session_number),
                data = group, 
                subset = (party_id != 4) & take_followup == 1) 
mfn.oft <- lm(other_ft ~ as.factor(treatment) + as.factor(session_number),
              data = group, 
              subset = (party_id != 4) & take_followup == 1) 
mfn.lt <- lm(less10 ~ as.factor(treatment) + as.factor(session_number),
             data = group, 
             subset = (party_id != 4) & take_followup == 1) 
mfn.gtf <- lm(gtfifty ~ as.factor(treatment) + as.factor(session_number),
              data = group, 
              subset = (party_id != 4) & take_followup == 1) 
mfn.ocft <- lm(other_cand ~ as.factor(treatment) + as.factor(session_number),
               data = group, 
               subset = (party_id != 4) & take_followup == 1) 
mfn.trait <- lm(avg_trait ~ as.factor(treatment) + as.factor(session_number),
                data = group, 
                subset = (party_id != 4) & take_followup == 1) 
mfn.trust <- lm(trust_other ~ as.factor(treatment) + as.factor(session_number),
                data = group, 
                subset = (party_id != 4) & take_followup == 1) 
mfn.sd <- lm(soc_dist ~ as.factor(treatment) + as.factor(session_number),
             data = group, 
             subset = (party_id != 4) & take_followup == 1) 

## table replicated on just the follow-up sub-sample 
stargazer(mfn.index, mfn.oft, mfn.lt, mfn.gtf, mfn.ocft, mfn.trait, mfn.trust, mfn.sd,
          se = starprep(mf.index, mf.oft, mf.lt, mf.gtf,mf.ocft, mf.trait, mf.trust, mf.sd),
          column.labels = c("AP Index","Out-Party FT","Out-Party FT< 10","Out-Party FT > 50","Out-Party Cand. FT",
                            "Out-Party Traits","Out-Party Trust","Social Distance (Reversed)"),
          covariate.labels = c("Homogeneous Discussion","Heterogenenous Discussion"),
          digits = 2, 
          keep.stat = c("n","rsq"),
          type = "html",
          out = "tables/appendix_effects_on_just_on_followup_sample.html")
